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By means of an adapted mean-field expansion for large fillings n S> 1, we study the evolution of 
quantum fluctuations in the time-dependent Bose-Hubbard model, starting in the superfluid state 
and approaching the Mott phase by decreasing the tunneling rate or increasing the interaction 
strength in time. For experimentally relevant cases, we derive analytical results for the temporal 
behavior of the number and phase fluctuations, respectively. This allows us to calculate the growth 
of the quantum depletion and the decay of off-diagonal long-range order. We estimate the conditions 
for the observability of the time dependence in the correlation functions in the experimental setups 
with external trapping present. Finally, we discuss the analogy to quantum effects in the early 
universe during the inflationary epoch. 
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I. INTRODUCTION 

The rapidly growing interest in the exploration of the 
dynamics of quantum phase transitions [l|| fosters our un- 
derstanding of the complex behavior of many-body sys- 
tems far from equilibrium. This concerns, in particular, 
the interplay of the microscopic degrees of freedom, their 
entanglement, and the resulting emergent behavior. The 
controlled study of the time development of correlation 
functions expressing these entanglement properties, and 
the dynamical emergence of correlations from an initially 
uncorrelated state has been undertaken for various sys- 
tems. The Bose-Hubbard model, describing the essential 
archetype of an emergence of strong correlations when 
one crosses a quantum phase transition point, was origi- 
nally introduced in a conventional condensed matter con- 
text, to explain certain properties of bosons in periodic 
and/or random potentials Its implementation with 
ultracold atoms [3| and the subsequent experimental re- 
alization of the superfiuid-Mott transition Q has caused 
a flurry of research activity. 

This activity is reviewed from a theoretical point of 
view in Q , while a number of recent experimental efforts 
studyingultracold atoms in optical lattices are covered, 
e.g., in [6|. Initially, theoretical studies were primarily 
devoted to the transition from the Mott regime, where 
number fluctuations are frozen (for commensurate filling 
of the lattice sites), to the superfluid side of the transi- 
tion for which, conversely, phase fluctuations are frozen 
0, S, B [13, [lH ■ Most of these investigations were done 
numerically, with the exception of certain exactly solv- 
able cases like the Ising chain in a transverse fleld 

We discuss here in detail a mean-field approach, first 
presented in which enables the (in some particu- 

lar cases analytical) rigorously controlled calculation of 
quantum correlations developing in rapid quenches from 
the superfluid to the Mott insulating phase. Such a 



number-conserving and hence controlled mean-fleld ap- 
proach is valid at large filling n 1 of the lattice sites, 
with the square root of the inverse fiUing l/^/n providing 
the expansion parameter Here, we supply in par- 

ticular analytical estimates for the applicability of the 
results obtained in [13!] to the experimentally relevant 
harmonically trapped case, by comparing the decay time 
of superfluid coherence to the propagation of the distur- 
bances in the system induced by the quench. It should be 
noted that the mean-field "hydrodynamic" limit of large 
site fillings considered in the following is not of purely 
academic interest. While many experiments on the Mott 
transition are carried out at small filling of order unity, 
experiments on number squee zing at large filling have in- 
deed been performed as well [iSl. [l^. [l7|. already on an 
early stage of research into the possible occurrence of the 
Mott insulator transition llal. 



II. THE BOSE-HUBBARD MODEL AT LARGE 
FILLING 

The Bose-Hubbard model describes interacting bosons 
on a lattice, hopping from site to site. Restricting our- 
selves to (two-body) contact interactions, i.e., that two 
bosons interact only if at the same site a, and to one 
single-particle state at each site, adding a one-particle 
scalar potential term, the Hamiltonian reads 



Q/3 



The coupling constant U is linear in the bulk contact in- 
teraction coupling constant g. The externally imposed 
scalar potential to additionally confine the atoms in the 
lattice is in most experimental situations to date to a 
goo d approximation harmonic, Va oc a^, or linear like in 
|17l |. Va (X a. The matrix Map describes the fact that, in 
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the lattice, the effective mass (and possible higher terms 
in a gradient expansion), can in general depend on posi- 
tion and the direction of hopping of the particles from site 
to site. For the simplest example of a one-dimensional 
chain with nearest-neighbor hopping, we have Map — 
5a,p - \ {5a,f}+i + 5a,i3-i + 5a-i,i3 + 5a+i,i3)-, and the effec- 
tive mass of the bosons is given by 1/m* — Ja^, where a 
is the lattice spacing. 

At large (average) filling n ^ 1 the Bose-Hubbard 
Hamiltonian |(T]) can be mapped to the so-called quantum 
rotor model (cf., e.g., |il!,lld|). Insertion of the (quantum) 
Madelung transformation da = e**^" V^hx into |(T]) and ex- 
pansion into inverse powers of n yields the quantum rotor 
Hamiltonian 




I J cos[(?!)a - 4>a+i] + -^{na- nY 



lJC0s[(l>a 



U 



(2) 



where we have taken for simplicity a ID lattice with no 
external trapping, Va = 0. As an important ingredient, 
we use in this representation that local number fiuctua- 
tion Sfia = ha — n and phase variables are conjugate in 
the limit of large n. 
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(3) 



The fact that the number and phase variables are canon- 
ically conjugate variables (implying the very existence of 
a phase operator) depends on applying the mean-field 
(and effectively hydrodynamic, i.e. coarse-grained) limit 
71 ^ oo has been known for a long time cf., e.g., J!). '2{\\. 
The problems arising without the limit n — > oo can be 
seen by means of the full commutator 



Pa,np 
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Since hp possesses a discrete spectrum with proper eigen- 
vectors I^q), taking the jnc) -expectation value of the 
above relation yields for a = /? a contradiction. In the 
large n-limit, however, the normalized fiuctuations Sha/n 
have quasi-continuous spectra and in this sense Eq. ([3]) 
provides a valid effective description. 

The form ((2]) provides a nice and intuitive understand- 
ing of the two phases: In the superfiuid state, the first 
term on the right-hand side dominates and (^q — 4>a+i is 
small (phase coherence). In the Mott phase, on the other 
hand, the second term wins and the number fluctuations 
Sha become small. Therefore, we may directly read off 
the scaling of the critical point Jc, which separates the 
two regimes Jc = 0{U/n). We depict a schematic repre- 
sentation of the superfluid-Mott transition on a quadratic 
2D lattice in Fig.[Tl 




FIG. 1: [Color online] Localization (superfluid-Mott) tran- 
sition of particles in an optical lattice. For increasing lat- 
tice depth and commensurate fllling (integer number of par- 
ticles per site), the system performs a quantum phase tran- 
sition from a superfluid declocalized phase with highly mo- 
bile bosons (left) to a localized insulator phase with immobile 
bosons at commensurate filling (right). The present analysis 
investigates the temporal evolution of correlations when one 
is rapidly going from the left (superfluid) to the right (Mott) 
side. For the sake of this representation, we display the non- 
mean-fleld case of unit fllling, n = 1, here. 



III. DYNAMICAL MEAN-FIELD THEORY 

We first consider the case without external trapping, 
Va = 0, for which the filling is homogeneous (ha) = n. 
From id]), the Heisenberg equations of motion for the 
lattice field operator are (we take h=l throughout) 



idtCla = J MapClp + Uhada 



(5) 



To proceed, we now employ a lattice version of the 
number-conserving mean-field expansion which in its 
continuum version has been introduced in (2ll . ^2f\ . Since 
the Hamiltonian cannot be diagonalized exactly, a 
controlled analytical approach requires a small or large 
parameter such that we can employ an expansion into 
powers of this control parameter. One option would be 
the weak-coupling regime J » C/, where an expansion 
into powers of U j J justifies the mean-field approach. 
However, since we want to study the sweep to the Mott 
phase, this ratio does not remain small and hence cannot 
be used as a control parameter. Therefore, we focus on 
the case of large fillings n 1 here and use l/n as a small 
expansion parameter (23l |. I.e., in the limit of large fill- 
ings n, the full site field operator is expanded into terms 
of various power in n. 



A 



(6) 



The operator A ~ a^.^a'^MY^^^I'^N^I'^ with ds = Yl,a 
accounts for the conservation of the total number of parti- 
cles N — A^A — J2a "a. Due to the number-conserving 
nature of the ansatz, we have (da) — exactly, as op- 
posed to non-number-conserving approaches, in which 
{da) = 0{y/n). We stress that a number-conserving 
formalism is essential to describe correlation functions 
accurately, especially those of higher order. 
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As stated before, the above mean-field expansion re- 
quires tiiat, at eadi site, the integral filling Ua = {ha) = 
{al^aa) e N is much larger than unity, 1. The idea 

of ([6]) is to expand the original operator da into pow- 
ers of no, ^ 1, considering the (formal) limit in which 
n t oo but with the chemical potential ^ = Un remain- 
ing finite and fixed, such that U = 0{l/n). The leading 
term in the above expansion ([6]) is the order parame- 
ter tpa — C'(-yn) describing the condensate part. The 
linearized quantum corrections Xa — 0{nP) are decom- 
posed into single-particle contributions, and correspond 
to Bogoliubov quasiparticle excitations above the super- 
fluid ground state, after a Bogoliubov transformation to 
a quasiparticle basis. For the vahdity of the expansion 
([6]), the remaining (non- linear) higher-order corrections 
Ca containing multi-particle contributions must remain 
small, Ca ^ 1, during the whole temporal evolution of the 
system. Since we are dealing with formally unbounded 
operators, the condition Ca 1 should be understood as 
a weak hmit, i.e., (Ca) <C 1 and (ClCa) ^ 1 etc. 

Thus, Cq, of order unity signifies that one approaches 
strongly correlated states, i.e., the Mott regime. These 
"strong correlations" are due to interactions between 
the quasi-particle excitations Xa , which become stronger 
when approaching the transition, and eventually gener- 
ate a gap due to non-perturbative effects. The emer- 
gence of this gap denotes the transition from the delo- 
calized atoms (superfluid state) to the localized atoms 
(Mott state), see, e.g., fH]. In what follows, we are in- 
terested in describing the approach to the localized state 
of the atoms, and how in a dynamical quantum phase 
transition the corresponding decrease in the phase order- 
ing takes place. 

The Bogohubov-de Gennes equations for the quasipar- 
ticle excitations at site a, which are by deflnition de- 
scribed by the first order operators Xa, read 

idtXa = Afa;3X/3 + 2U\i>l\xo. + U^lxi , (7) 



with Vo being determined by the solution of the Gross- 
Pitaevskii mean-field equation 

zStV'o = Afa/3?Ao + C/|i^o|Vo = C/IV'olVo , (8) 

where the final equality holds for a spatially homogeneous 
ground state. All the residual terms, remaining after the 
insertion of ^ into |[5|), determine the higher-order cor- 
rections Ca to the mean-field expansion which then 
evolve according to the equation of motion 

idtCa = J^Ma/3C,3+2;7|V^^|^^^^^^2^t + (9) 

13 

+2UiJoxixa + Ur.xl + Uxixl + 0{UU) . 

Deep in the superfiuid phase (our initial state), the 
higher-order corrections Ca are small and the mean-field 



expansion ^ works very well. If we approach the Mott 
phase, however, these corrections start to grow accord- 
ing to Eq. ^ and, at some point, the mean-field ex- 
pansion ([6]) breaks down. The characteristic time-scale 
of this breakdown can be estimated from the nonlin- 
ear source terms in Eq. ([9]), which are suppressed to 
0{\/y/n) in view of the imposed constancy of the chem- 
ical potential, so that U = 0{l/n), ipg = 0{^/n), 
and Xa = 0(nP). Therefore (starting in the superfiuid 
phase) , the higher-order corrections remain small as long 
as we have Ut^/n 1, i.e., for evolution times which 
are of order t — 0{yjn). Thus, while Bogoliubov mean- 
field theory is, in principle, also valid (initially, i.e., for 
J ^ U) for small values of n = 0(1), the Ca corrections 
in ^ , which describe correlations beyond those derivable 
from the mean-field plus single-particle fiuctuations, i.e., 
from -00 and Xa in ([6]) , begin to grow very quickly in off- 
equilibrium situations, then invalidating the Bogoliubov 
approach at small filling. For large fillings n ^ 1, how- 
ever, the mean-field expansion ([6]) can be extrapolated 
to relatively long time-scales t = 0{^/n), which allows 
us to study the sweep analytically (with Ca serving as an 
indicator for the vahdity of the mean-field approach) . 

Linearizing the polar decomposition of the fundamen- 
tal field operator Oq, = exp{i(/)Q}\/n^, we may identify 
the fiuctuations Xa — iJolSha/i^n) + iS(j)a] + ©(l/y^) 
in terms of the linearized number fiuctuations Sha and 
the conjugate phase fiuctuations (5(/)q according to |[3|. 
Eq. ^ can be diagonalized by a normal-mode expansion 
into the eigenvectors Mapv^ = Xk,v" of the hopping ma- 
trix Map, with the eigenvalues labeled by the general- 
ized momenta k. To this end, we expand the number and 
phase fiuctuations via Jn^ — v"5ha and 54)k = v"54)a, 
respectively, and insert them into 

Xc. = V^e-'''^'^''\6ha/2n + i6^a), (10) 

where we have factored out the time-dependence of 
the Gross-Pitaevski! wave function with the appropriate 
phase factor V'o = \fne^"'' ^ as follows from Thus, 
we obtain from the Bogoliubov-de Gennes equations ((?]) 
two real equations for the time evolution of density and 
phase fiuctuations in the eigenbasis labeled by k, 

dM. = - + [/) a. , (11) 

for generally time-dependent J and U . In case that both 
J and U are constant in time, the above equations result 
in the well-known Bogoliubov spectrum [H, HB] 

^ j2x2 ^ 2nUJK ■ (12) 

We now discuss two possible routes to approach or 
cross the phase transition from the superfiuid to the Mott 
side dynamically. Either one decreases the tunnehng rate 
in time or, alternatively, the interaction is increased to 
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suppress the superfluid density at given filling and thus 
cross the transition line. It is demonstrated that in sev- 
eral particular cases for J — J{t) or U — U{t), respec- 
tively, the Bogoliubov-de Gennes equations ([7]) can be 
solved analytically. 

A. Decreasing the tunneling rate 

Combining the two equations (fTTj) , we get the following 
second-order equation for the number fluctuations [anal- 
ogously for phase fluctuations, see Eq. l(26|) below] 

/did \ 

f ^ J ^ + A« [ JA, + 2Un] \ Sn, = . (13) 

Dividing this equation by J and A^, and deflning a new 
time coordinate depending on the mode index k, 

dT^ = \^Jdt, (14) 

we obtain an equation containing the operator d'^/dr^, 

/ 92 U2n\ ^ 

We note that this equation now contains the ratio U/ J 
only, but not U and J separately, cf. Sec. IIIICl below. 



Let us first study the case J = J{t) while U = const. 
Even though it is not possible to give a closed solution 
of Eq. (fTSl for arbitrary time-dependences J{t), there 
are several cases which do admit analytic expressions in 
terms of Hankel H^}^ or Whittaker Wj/,^ functions [3l] |. 
We Hsted a few cases in table H] below. 

In view of the asymptotic M_oo behaviour of the Han- 
kel and Whittaker functions [3l|, we see that the number 
fluctuations Shu oscillate forever in the last three cases 
{J{t) oc ,t~^^^,t~^/^) though with a decreasing am- 
plitude and frequency. In the case of an exponential - i.e., 
much faster - sweep J{t) = Joe"^*, on the other hand, 
the solutions Suk, do not have enough time to adjust to 
the externally imposed change of J{t) and freeze at a fi- 
nite value (non-adiabatic behaviour). Finally, the second 
case J{t) = Jo{jt)~'^ just marks the border between the 
two regimes (eternal oscillation versus freezing). Con- 
sequently, the number fiuctuations Sh,^ vanish for late 
times in this situation (for Ak > 0). 

The asymptotic behaviour can be interpreted nicely in 
terms of the analogy to cosmology sketched in Sec. IVl 
below. The freezing of the number fiuctuations (Sn^ 
for J{t) = Joe~^* can then be understood via the 
emergence of a horizon analogue - whereas for J{t) oc 
t"^/"^, such a horizon is absent. The critical 

behaviour J{t) — Jo(7t)~^ precisely marks the limit for 
horizon formation, cf. Sec. |Vl 



Jit) 


Solution 


Argument 


Indices and constants 


Asymptotics t 1 00 


Joe--" 




Xk = Tk = —Jo)^K.e~'^^/j 


H = 1/2, ly = [/n/7 


xn T 






Xk = = -JoAK(7t)~V7 




T 




hI'\x^) 


x^ = 2V2(7nJoA«(7t)i/V7 


V = ±2jJoA„/7 






XK^^^Wi^^^iix^) 


x^ = 3V2(7nJoA.(7t)^/V7 


= 1/16, = 9(JoA,)V(32C/n72) 






^hI'\c^xT) 


x^ = Vt* + JoA«/(2C/n) 


u = 1/3, c„ = V2l7nJoA./(37) 


x^^ 00 



TABLE I: A few analytically solvable cases for time-dependent J — J{t). The last column indicates the behavior of the 
argument Xk at late laboratory time. 

I 

very quickly. For the Bose-Hubbard model, this time- 
dependence would allow us to approach the Mott state 
very efficiently. 

However, in an experimental realization, a dynamics 
J{t) — Jo{'jt)^'^ requires some fine-tuning and is proba- 
bly hard to achieve. The exponential sweep is more inter- 
esting from a theory point of view (since it yields non-zero 
frozen number fluctuations) and should also be a better 
approximation to a realistic experimental situation. The 



It is very illustrative to compare the various cases 
discussed above to a harmonic oscillator with a time- 
dependent damping term and/or spring constant. In the 
last three cases {J{t) cx t^^, t^^/'^, t~^/^), the spring con- 
stant dominates (underdamped oscillator) while the ex- 
ponential sweep J{t) — Joe^''* induces a transition to 
the overdamped regime at some time. The boundary 
case J{t) = Jo(7i)~^ would then correspond to criti- 
cal damping, where the solution Suk approaches zero 
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experimental relevance becomes apparent considering the 
fact that, in sufRciently deep d-dimensional simple cu - 
bic lattices, J{t) cc {Voit)/ER)^/^ exp{^2^Vo(t)/E^} 
and U oc {Vo/Ei^y^^ hold, where Voit) is the time- 
dependent lattice depth given by the laser intensity and 
Ejf = 7r^/(2ma^) is the (constant) recoil energy, with 
TO the bare boson mass |33|- Apart from logarithmically 
slow corrections, an exponential sweep J(t) oc e~'^* there- 
fore corresponds to increasing the laser amplitude (and 
therefore ^/Vq) linearly in time, with U then remaining 
approximately constant. 

For the case of an exponentially decreasing tunneling 
rate, a universal "scaling" solution exists: As we may 
infer from the table HI the solution (i.e., the indices fi and 
V of the Whittaker function W^,^) then do not depend 
on At. The "scaling time" from l|14p . depending on the 
mode index k, reads t„ — ~Xi^J{t)/'y, so that Eq. (fTSl is 
transformed into a scale invariant equation of the form 
of Eq.dlll), 



1 



2 Un 



0. 



(16) 



The only remaining dimensionless parameter determin- 
ing the relevant universality class of solutions of this 
equation \?,v = Un/^. This parameter, equal to the ratio 
of chemical potential ^ = Un (i.e., internal energy scale) 
and sweep rate 7 (i.e., external time scale), represents a 
measure of the rapidity of the externally imposed sweep: 
V ^ 1 implies a slow and ^ 1 a fast (non-adiabatic) 
sweep. The adiabaticity parameter v determines the na- 
ture of the final state. Starting in the superfiuid phase 
and ramping down the tunneling rate very rapidly, <C 1 , 
the system will have no time to react to this change. For 
later times, the vanishing hopping J = prevents an 
equilibration of the number fluctuations, i.e., they will 
be as large as in the initial coherent state in this situa- 
tion. The slower we sweep J(i), the more time the system 
has to adapt to this external change and the closer the 
final state will be to the Mott state (i.e., smaller number 
fluctuations will result from a slow sweep). 

As listed in table HI the analytical solution of flGl) can 
be found in terms of the Whittaker functions Wi^,ii2 
It allows us to obtain the following exact Bogoli- 
ubov transformation from bare density /filling fiuctuation 
to initial vacuum quasiparticle operators 6^, labeled by 
the mode index k. 



5n^ = V^e--''/^ W^,,a/2(2ir,)feK 



(17) 



The initial vacuum quasiparticle operators annihilate 
the adiabatic superfiuid ground state 6^ I in) = at early 
times Tk [ —00, where the modes oscillate like e^*"^". 
(Note that the scaling time t„ is dimensionless). The 
phase fiuctuations are obtained using the relation (50k = 
— ^dSn^/dr^ following from ifTTj) and (fT4| . 



Ok 



dr.. 



h.c. 



(18) 



The analytical scaling solution thus obtained represents, 
to the best of our knowledge, the first example of an exact 
solution of the (nonintegrable) Bose-Hubbard model on 
the Bogoliubov mean-field level in a dynamical situation. 

Due to the perfect scaling solution in Eq. lfT7|) . the 
frozen value of the number fiuctuations at late times, 
when — *■ 0, is independent of k, but the decaying cor- 
rections do depend on the eigenvalue A„: 



(Shi) = (in| ((5n„)^ |in) = n 



1 - e- 



2'KV 



0(tAKe-^*)(19) 



Since the leading term is independent of k, it just yields a 
local (oc 5a, 13) contribution after the mode sum (k a) 
and thus leads to frozen on-site number variations. 

In contrast to the number fiuctuations which freeze 
according to ifTQ]). the conjugate phase fiuctuations grow 
(as one would expect when approaching the Mott phase) . 
From our analytical result (fT8|) . we conclude that they 
increase (initially) quadratically in time: 



1 -e- 



2'Kn 



(20) 



Again, like for the number fiuctuations, the leading (first) 
term is independent of the mode index k and yields the 
on-site phase fiuctuations {Scf)^)- The off-site phase cor- 
relations {54>aS4>f3) ^ corresponding to the second term, 
grow linearly in time (initially). 



B. Final state 

All the results so far were obtained by a controlled 
extrapolation of the mean-field expansion ^ from the 
weak-coupling (J dominates) into the strong-coupling 
regime {U dominates) and hence are only valid as long 
as the quantum depletion (xlxa) is small, i.e., the con- 
densate ipo dominates. In the strong-coupling regime, 
however, the quantum depletion grows (on a time-scale 
of order ^Jn) and finally invalidates the mean-field expan- 
sion ([6]). Fortunately, we may also analytically describe 
the ensuing stages of the quantum evolution, because the 
tunneling rate J{t 3> I/7) ^ 1 is exponentially small 
and can be completely neglected. In this completely 
interaction-dominated limit, the evolution of the site op- 
erators can be approximated by dda/dt = —iUfiada, 
which possesses the simple exponential solution 



ia{t) — eyip{~iUn^t}c 



(J « 1) . 



(21) 



Consequently, we may calculate the time evolution of 
correlation functions throughout the dynamics (i.e., for 
all t) by using the results of the mean-field expansion, 
which are valid for intermediate times with ^ 1 and 
Uty/n ^ 1, as initial conditions, and then switching to 
the above solution (|2T|) for later stages. 

A frequently used indicator for distinguishing the su- 
perfiuid from the Mott phase is the off-diagonal long- 
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range order (ODLRO) usually associated to the first- 
order correlation function {al^{t)a/3{t)) . Using the re- 
sults above, we may derive an analytical expression for 
the time-dependence of (ajj(t)a^(t)) during and after the 
sweep. From (|2T|) . we obtain for the correlator 



exp{iU{na-ni3)t}) +OiV7i) . (22) 



On the other hand, the frozen first-order number fiuctu- 
ations Sfia are then in a squeezed state which can (for 
n 1) be approximated by a (continuous) Gaussian dis- 
tribution. For a Gaussian variable X with {X) — 0, the 
exponential average yields (exp{iX}) = exp{— (X^)/2} 
and hence we finally get 



{dl{t)af){t)) « nexp{^ 



(23) 



where the on-site number variations are given by ([19 



A\n^) = (nl) - {n^f = (Shi) 



2tiv 



(24) 



The result l|23p represents an extension to the present 
period lattice geometry of the result obtained in [23 |. 
where the interaction-induced decay of coherence in a 
double-well trap was studied. The Fourier transform 
of the first-order correlation function {al^{t)ap{t)) deter- 
mines the corresponding momentum distribution func- 
tion gi[k). The decay of the off-diagonal long-range order 
(ODLRO) in (|23l) thus directly corresponds to a tempo- 
ral decrease of the peak in pi(k) at fc = 0, measurable in 
time-of-fiight experiments |29l. l30l|. 

In summary, the state obtained with J{t ^ I/7) ^ 1 
(while still maintaining Ut^/n <C 1, i.e. at intermediate 
times), which has frozen number fluctuations according 
to (fT9|) , forms an appropriate initial many-body quantum 
state for the further evolution beyond mean-field into the 
strongly correlated Mott phase. 



C. Increasing the interaction coupling 

Another possibility to approach the Mott phase dy- 
namically is to increase U in time. This can be experi- 
mentally reaHzed using a time-dependent sweep through 
a Feshbach resonance, which varies the s-wave scattering 
length as, and thus U only, while keeping Vq and thus 
J fixed. Since Eq. (fT5|) depends on the ratio U/ J only, 
increasing U is analogous to decreasing J. Therefore, we 
can establish an exact "duality" between the results ob- 
tained for J = J{t) and the present U = U{t), i.e., every 
analytic solution in table U of Sec. IIII Al corresponds to 
a dual expression for U = U{t) after incorporating the 
transformation t ^ t oi the time coordinate in Eq. (fT4l) . 
For a power-law dependence U on t"' , Eq. I|15p possesses 
formally the same solution as for J (x t~". Transforming 
back to the laboratory time via Eq. ([14]), this corresponds 
to J cx while the two limiting cases a = —1 

and a = 00 correspond to exponential dependences, see 



At) 


Dual Uit) 


Dual Range 


Joe-'" 




-00 < t < 






-00 < t < 






—00 < t < 00 






< t < 00 




Uo-ft 


< t < 00 



TABLE IL Important dual cases for temporal variations of 
U and J. Given the solution for the first column, we can 
immediately conclude, by a simple time transformation, on 
the solution for the second column, and vice versa. 



Tab.lTTl Due to the transformation t r of the time co- 
ordinate in Eq. I|14p , the solutions for U (t) do not freeze 
in terms of the laboratory time t ^ 00. 

Again, it is likely that most of the dynamics U {t) are 
hard to realize experimentally. Therefore, we focus on 
the linear increase, U{t) = jt, in the following, since this 
case is probably close to an experimental set-up. For a 
linear growth (dually corresponding to J{t) oc 1/^/t), we 
can define a scaling time using a simple K-dependent shift 
of the time origin. 



= t + - — dr^ = dt , 

2nj 



(25) 



Introducing this scaling time into the equation for phase 
fluctuations 



d 



d 



JX 



dt 2nU{t) + J\^ dt 
we obtain a scaling equation of the following form 
d 1 d 



, (26) 



+ 2^7 JA« U0« = 



(27) 



The solution of this equation can be found in terms of 
Hankel functions and leads us to the following Bogoliubov 
transformation for the phase fluctuations 



-h.c. (28) 



where Hi''^ and 



{hI^'')* are Hankel functions of 
the flrst and second kind, respectively. Note that the 
index is mode number k- independent, like the solution 
for exponential decrease of J in l[T7j) . In contrast to the 
above-discussed case of an exponential decay of J, ap- 
proaching the "hard-core" limit by letting U grow (lin- 
early), we do not obtain frozen number fluctuations at 
late times ^ 00. We then have, instead, increasingly 
rapid oscillations of fluctuating flUing 



1 d5(t>^, 



1/2 



h.c, (29) 
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where the operators (7^ are proportional to the in 
(|28l) . Due to the asymptotic (large x) behavior H^{x) 
y/2/nx exp[i(x — — 7r/4)], we have at late times 

Sn^ cx exp[i| ^/2n'YJX^t^^'^]+h.c. From ([28]), we 

derive the corresponding slow increase of the amplitude 
of phase fluctuations, (5(/)„ cx 

h.c, with the same increase of the oscillation frequency. 



D. Quantum depletion 

Approaching the Mott phase entails that the conden- 
sate becomes depleted due to Bogoliubov quasiparticle 
excitations created above the (still superfluid) ground 
state. We stress that the quantum depletion we calculate 
below explicitly depends on J = J{t) [or U = U{t)] in a 
particular manner. Thus it is possible, within our con- 
trolled mean-fleld scheme, to unambiguously identify the 
number squeezing (related to the depletion) caused by 
the creation of the excitations, and to distinguish it from 
effects coming from the time dependence of the mean- 
fleld in the Gross-Pitaevskii equation (which were dis- 
cussed in [38]). 

The (relative) depletion is defined by the expression 



The continuum scaling equations of motion in ID read 



V 



(30) 



Concentrating here on J = Joexp[— 7^], where as shown 
in section lTlI AI number fluctuations freeze, and using that 
the local depletion is given by XaXa = Sn^/ {4:n)+n5(j)'^ + 
1/2, we have from ifB ]) and 



V 




> 1 



1 

2n 
(31) 



which is valid provided that the mean-field condition 
Uty/n <C 1 holds - which implies P <C 1. The phase fiuc- 
tuations caused by the sweep always dominate the num- 
ber fluctuations in the limit assumed throughout, n ^ 1 
(while still keeping Ut^/n <C 1). 



E. Comparison to continuum scaling 

Up to now we have treated the homogeneous case. If 
the lattice is embedded in an external harmonic trap- 
ping potential, like in experiment, the derived analyti- 
cal solutions for J cx exp[— 7t] and U (x t are not ex- 
actly valid anymore. However, in the continuum limit, 
the well-known scaling approach for a bulk gas in a har- 
monic trapping potential [33, 34, 35] can be applied and 
it is rather interesting to contrast this approach with our 
scaling solutions on the lattice ifTTj) and (|29l) . 



bit) 



■m 

Jo 



m 
j{t) 



bit) 



Jit) Uit) 



Jo Uo b-^HY 



(32) 



where bjt) denotes the time-dependent scaling parame- 
ter jsl, [3, [3^ while cjQ, Jo, and Uo are the initial values 
for trap frequency, tunneling, and interaction rate, re- 
spectively. We conclude from Eq. (|32)) that the temporal 
change of Jit) corresponds to a change in the (inverse) 
effective mass of the bosons. 

For Jit) = Joe~'^*, the effective harmonic trapping fre- 
quency becomes exponentially slower, and the scaling fac- 
tor motion consequently overdamped due to the resulting 
constant (Ohmic) damping. This damping corresponds 
to the derived freezing of the number fluctuations and 
increase of phase fluctuations on the lattice: 



bit) + 76(t) + e-^'ulbit) = 



b^it) 



(33) 



Conversely, for U it) linearly increasing in time, and con- 
stant J, due to the driving term on the right-hand side 
of l(32|) , the scaHng factor begins to oscillate increasingly 
fast, like already observed from the exact solution on the 
lattice, Eq. l(29l) . 



IV. 



NONEQUILIBRIUM QUANTUM EFFECTS 
IN INHOMOGENEOUS SYSTEMS 



In the additional presence of an optical lattice, the be- 
havior of the system depends on the relation between 
the inhomogeneity of the trap potential Va and the cen- 
tral filling ncontcr ^1. If the potential Va is rather 
shallow, the system develops a "wedding cake" structure 
near the boundaries, where the filling is small and hence 
the Mott phase emerges [1, [3H[. For stronger inhomo- 
geneities of Va and/or larger central fiUings ncontcr ^ 1, 
on the other hand, two-body interactions C/(aJj,)^a^ in 
(HJ) can be neglected in comparison with the tunneling 
term JMc^aJ^a^ and the potential gradient of Vafia due 
to U = ©(l/nccntcr)- Intuitively speaking, the various 
rims of the "wedding cake" would be compressed into a 
single lattice site and hence disappear. 

Lowering J towards the (first) phase transition point, 
the wedding cake structure gradually develops by Mott 
insulator shells propagating inwards from the outer low 
density edges of the gas cloud. Therefore, we have to 
study the question of whether the results derived above 
for the homogeneous case - such as the temporal decay of 
the coherence peak of the momentum distribution func- 
tion gi(fc) due to number squeezing effects - can still be 
observed in the presence of an external trap. Clearly, 
near the boundaries inhomogeneity effects will be impor- 
tant. However, in the central region, the homogeneity 
assumption should provide a good approximation - pro- 
vided that the momentum distribution function gi(fc) de- 
cays fast enough, i.e., before the Mott insulator shells 
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propagating inwards reach the center. From Eq. I|23p . we 
conclude that the typical time for the peak at zero wave 
vector of the momentum distribution function gi{k = 0) 
to decay (say, one e-folding) is generally scaling with 
1/y/n, i.e., we have 

U^-^ (.«1), t,^^^ (.»1).(34) 

To clarify the question above without resorting to brute 
force numerical computation, using simple analytical 
means, we consider in what follows various possible 
modes of the disturbance propagation of the corrections 
to mean-field, to determine the dominant, i.e. most 
rapid, mode of propagation. 



It is a well-established fact that shock waves in a homo- 
geneous medium can propagate much faster than sound; 
we therefore need to estimate if shock waves might over- 
take the sound waves in the proceeding quench. However, 
according to the study of [39| , the propagation of shock 
waves in the optical lattice proceeds at speeds slower than 
that of the sound waves for all values of J < the dis- 

persion relation l(35|) then has negative curvature for all k. 
This condition translates into {J/ Jc)(l/n'^) < ^. There- 
fore, close to the Mott transition, and for sufficiently 
large filling, "shock" waves always propagate slower than 
sound, as opposed to the uniform system. 

B. Propagation of the phase boundary between 
superfluid and Mott phases 



A. Propagation of sound and shock waves 

The speed of sound in a ID optical lattice, for large 
fining n, can be derived from the dispersion law js^l 



uil = 4nC/ J sin^ 



4 J2 sin'' 



ka 



(35) 



where the lattice spacing a = A/2 is determined by the 
laser wavelength A [the eigenvalues of Mafs in this case 
are Xk — 2sin^(fca/2), from which, using (fT2|) . the above 
dispersion law follows]. This leads for /c — > to the sound 
speed 



W nil J 



nUJ 
2mER 



(36) 



A condition for the undisturbed observability of the de- 
cay of gi{k) is, then, that the time for sound pertur- 
bations to propagate to the center is much larger than 
the decay time td, i.e., RtfI{cs)o ^ ^d, where the 
Thomas- Fermi size o f the harmonically trapped cloud 
Rtf = \/2mUnQ/Lo'^, with no the central fiUing and cen- 
tral sound speed (cs)o — Cs{J — Jo,n — uq). We then 
obtain that, using Jc^U/n 




TTUJ y j/Jc n 



has to be fulfilled for a rapid quench, v -^1. This condi- 
tion is practically always fulfilled if the optical lattice is 
not too tightly harmonically trapped: With a laser wave- 
length of A = 985 nm, we have En ~ 3.7 kHz for ^llb and 
Eft ~ 8.9 kHz for Na, while typical values for the on-site 
interaction are U ~ ■ ■ ■1Q~^ Er furthermore, 

Lu in the weakly confining direction(s) does not exceed 
100 Hz. For a slow quench, t he ri ght-hand side of the 
above inequality ((37| reads \fvjn (respectively ^Jv/n) 
instead, and the inequality will only be fulfilled if v is 
not too large, i.e., if the sweep is not too slow (as one 
would expect). 



In the two examples above (sound and shock waves), 
we studied the propagation within a given phase (the 
superfiuid phase). However, one might object that the 
motion of the superfiuid-Mott phase boundary could per- 
haps be much faster. Clearly, the typical time scale to 
overcome one lattice site will again be set by the tunnel- 
ing rate J, but this rate could possibly be enhanced (or 
suppressed) by a large filling factor n 3> 1. In general, the 
correct description of the propagation of a phase bound- 
ary in the Bose-Hubbard model is a very interesting and 
quite involved problem. In the following, we derive a 
rough estimate using the simpler two-site Bose-Hubbard 
model (a Bose-Einstein condensate Josephson junction), 
described by the Hamiltonian 



H = J{aia. 



a\a2) 



0{U). 



(38) 



In one site (ai), we model the superfiuid phase by a co- 
herent state, while the other side (02) is prepared in a 
number eigenstate simulating the Mott phase, cf. Fig.[2l 
Therefore, the combined initial state factorizes 



(39) 



into a coherent state and a number eigenstate jV'^) 

ai 10 = ai |^°) , n2 |^2°) = ^2 |^2°) ■ (40) 

Now the question is the following: how fast does the 
phase coherence of the coherent state (superfiuid) decay 
due to the coupling with the other site (Mott)? Since 
the propagating superfiuid-Mott phase boundary is (like 
a shock wave) presumably very sharp, this site-to-site 
behavior should provide a rough stimate for the prop- 
agation speed. Of course, there will also be an on-site 
dephasing caused by the self-interaction U, but this pro- 
cess precisely corresponds to the homogeneous growth of 
the phase fiuctuations considered in the present Article 
and is independent of the moving phase boundary. 

Without loss of generality, we assume ai € R, i.e., 
(pi)o = and (gi)o 7^ in terms of the canonical vari- 
ables qi = (a[ + di)/\/2 and pi = i{a\ — di)/\/2. The 
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^<q> 

A(p 




FIG. 2: [Color online] Phase space diagram (not to scale) of 
the coherent state (left) and the number eigenstate (right). 



initial phase uncertainty is given by, cf. Fig. [2] 



((A^i)'>o 



(9i>^ 







(4) 









, (41) 



which is the usual result for a coherent state. Now, if we 
insert the time evolution in the Heisenberg picture 

pi{t) - pi(0) + Jtg2(0) + 0{jh'') , (42) 

we find a gradual increase of the phase variance 



((A^i)')W = ((A^i)^)o + J^i^j% 



(43) 



which just corresponds to the dephasing induced by the 
coupling to the other site. Since the phase fluctuations 
of the other site are maximally large (gDo = 0{n2), we 
get together with (gi)o — 0{^/ni) the final result 



111 



(44) 



Hence, for equal (mean) fillings ni — n2 — n, the typical 
dephasing time scale ((A(pi)^) — 0(1) is just determined 
by the tunnehng rate J without suppression or enhance- 
ment by additional factors of n ^ 1. 

To summarize, the propagation velocity of the two-site 
model for the phase boundary is of order V2 = 0{aJ). 
The ratio of the maximal sound speed from (|36|) to V2 is 
given by 



V2 



(45) 



It follows that sound waves are the dominant, i.e., most 
rapid mode of propagation of the disturbances caused by 
the Mott insulator shells, due to the n-dependent collec- 
tive enhancement factor, which is absent at least from 
the simple two-site model for the phase boundary prop- 
agation. 



V. EFFECTIVE SPACETIME AT LARGE 
WAVELENGTHS 

The freezing of the number fluctuations obtained in 
section IIII Al can nicely be explained via the analogy to 



the kinematics of quantum fields in gravity, valid at large 
length scales and low energies. At large wavelengths 
A a, the lattice structure is not important for the 
propagation of disturbances and the system behaves like 
an ordinary (possibly inhomogeneous and moving) su- 
perfluid. [This is precisely the limit in which the scaHng 
Eq. l(32|) holds.] In the continuum Hmit, the excitations 
(phonons) possess a linear spectrum at low energies and 
behave in complete ana logy to scalar (quantum) fields in 
curved space-times [13, l4l| . 

The case of a decreasing tunneling rate J{t) is anal- 
ogous to an expanding universe. This can be seen by 
means of the evolution equation l|13p in the limit of small 
wavenumbers k — -v/SAT/a 



d 



3t 



0. 



(46) 



which is identical to the wave equation for a scalar field 
mode in an expanding universe [i^ . [i^ ]. If the speed 
of sound Cg(t) = Una?J{t) decreases fast enough (e.g., 
exponentially) , phonons cannot propagate arbitrarily far 
but may only travel a finite distance. In analogy to cos- 
mology, this corresponds to the emergence of a horizon 



Ah(t) = / dt' csit') =a dt' ^JJ{t')Un 



(47) 



The convergence of the above integral for t' — > oo then 
indicates the emergence of a horizon. This observation 
matches the conclusions of Section ITlI A| where we found 
that the solutions for J cx t~°' with a < 2 do not freeze 
at late times but oscillate forever (no horizon). For an 
exponential sweep, on the other hand, the integral con- 
verges and a horizon exists. Because its size shrinks expo- 
nentially Ah(i) = 2V JoUn exp[— 7i/2]/7, it engulfs any 
given mode with wavelength A after some time. There- 
fore, the evolution of the phonon modes passes through 
three stages: oscillation A ^ Ah(t<), horizon crossing 
A — Ah(i=), and freezing A ^ Ah(i>). Since each mode 
can be mapped to a harmonic oscillator, this evolution 
corresponds to the transition from underdamped to over- 
damped regime, which can be seen by rewriting l(46l) as 

(^^-j^^-J{t)Una'^^^5h{x,t) = Q. (48) 

Hence, identifying J / J with the Hubble constant (which 
is indeed a constant for the exponential sweep) , the freez- 
ing of the number fiuctuations and the creation of a num- 
ber squeezed state by the exponential sweep considered 
in this Article is completely analogous to cosmic inflation 
[12]. During this very early epoch of the cosmic evolu- 
tion, the rapid expansion of space induced a squeezing 
of the quantum fluctuations of the inflaton scalar fleld 
- traces of these frozen and amplifled quantum fluctua- 
tions can still be observed today in the anisotropics of 
the cosmic micro-wave background radiation. 
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VI. CONCLUSIONS 

For the Bose-Hubbard model with generally time- 
dependent coefficients J{t) and U{t), we developed a 
rigorously controlled and number-conserving mean-field 
expansion for large filling n ^ 1 as a generalization of 
the original theory of Bogoliubov to the lattice. This 
allows us to study non-equlibrium quantum phenom- 
ena occuring in the sweep from the superfiuid to the 
Mott phase. For two experimentally relevant cases - 
an exponential decay of the tunneling rate J{t) and a 
linear increase of the interaction coupling U (t) , respec- 
tively - we found exact scaling solutions which facilitate 
fully analytical expressions for the time-dependence of 
the Bogoliubov excitations and the resulting quantum 
depletion. Moreover, we were able to establish a du- 
ality between various power-law behaviors of U{t) and 
J{t), which enables transferring one solution obtained for 
J = J{t),U = const, into one for U = U{t),J = const., 
and vice versa. 

In the first case (exponentially decaying tunneling 
rate), we observe freezing and squeezing of the initial 
quantum number fiuctuations in complete analogy to cos- 
mology. The final state and its energetic position be- 



tween the (initial) coherent state and the Mott phase 
(final ground state) depends on a single adiabaticity pa- 
rameter ly, which is given by the ratio of the (external) 
sweep rate and the (internal) chemical potential. This 
parameter also governs the decay of off-diagonal long- 
range order. 

Since the Bose-Hubbard model is considered a proto- 
typical example for quantum criticality, we expect our 
findings to contribute to the general understanding of 
dynamical quantum phase transitions. Furthermore, by 
estimating the applicability of our results derived for the 
homogeneous case to the real experimental situation with 
harmonic trapping present, we found that the predicted 
effects should be observable with current optical lattice 
technology. 
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